#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include <limits.h>
#include <assert.h>
#include <mpi.h>

#include "halofinder.h"

int read_GADGET_Page(long PageNum, long *NumBuffParts, float **buffPx, float **buffPy, float **buffPz, float **buffVx, float **buffVy, float **buffVz)
{
  static FILE *fp;
  long k;
  char basefname[MAX_FILENAME];
  char fname[MAX_FILENAME];
  static float *partChunk,sqrta;
  static long chunk,TotNumPageParts;
  static long NumPartsInChunk,BaseNumPartsInChunk,NumChunks,CurrNumParts;
  static long InternalPageNum = -1;
  static fpos_t pos_p,pos_v;
  
  struct io_header_1 {
    unsigned int npart[6];
    double   mass[6];
    double   time;
    double   redshift;
    int      flag_sfr;
    int      flag_feedback;
    unsigned int      npartTotal[6];
    int      flag_cooling;
    int      num_files;
    double   BoxSize;
    double   Omega0;
    double   OmegaLambda;
    double   HubbleParam;
    char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  // fills to 256 Bytes
  } header1;
  
  int dummy;
  
#define SKIP fread(&dummy,sizeof(dummy),(size_t) 1,fp);
  
  if(PageNum != InternalPageNum)
    {
      InternalPageNum = PageNum;
      
      fp = fopen(haloFinderData.SnapshotFileList,"r");
      assert(fp != NULL);
      fgets(basefname,MAX_FILENAME,fp);
      fclose(fp);
      
      fp = fopen(basefname,"r");
      if(fp == NULL)
	{
	  sprintf(fname,"%s.%d",basefname,(int) PageNum);
	  fp = fopen(fname,"r");
	  assert(fp != NULL);
	}
      
      SKIP;
      fread(&header1,sizeof(header1),(size_t) 1,fp);
      SKIP;
      
      TotNumPageParts = 0;
      for(k=0;k<5;k++)
	TotNumPageParts += (long) (header1.npart[k]);
      
      sqrta = (float) (sqrt(header1.time));

#ifdef DEBUG
#if DEBUG_LEVEL > 1      
      fprintf(stderr,"NumPart = %ld \n",TotNumPageParts);
      fprintf(stderr,"header1.npart[1] = %u \n",header1.npart[1]);
#endif
#endif
      
      BaseNumPartsInChunk = (long) (IO_CHUNK_SIZE_MB*1024.0*1024.0/((double) (sizeof(float)))/3.0); //IO_CHUNK_SIZE_MB MB read chunks for parts
      partChunk = (float*)malloc(sizeof(float)*3*BaseNumPartsInChunk);
      assert(partChunk != NULL);
      NumChunks = TotNumPageParts/BaseNumPartsInChunk;
      if(NumChunks*BaseNumPartsInChunk != TotNumPageParts)
	NumChunks += 1;
      
      chunk = 0;
      CurrNumParts = 0;
      
      /*
	will keep track of position in file between calls of function using the fgetpos and fsetpos
	use one file position for positions and another for velocities
	need to update each position pointer after the read of the file
      */
      SKIP;
      fgetpos(fp,&pos_p);
      fseek(fp,(long) dummy,SEEK_CUR);
      SKIP;
      
      SKIP;
      fgetpos(fp,&pos_v);
    }
  
  if(CurrNumParts + BaseNumPartsInChunk > TotNumPageParts)
    NumPartsInChunk =  TotNumPageParts - CurrNumParts;
  else
    NumPartsInChunk = BaseNumPartsInChunk;
  
  //alloc mem
  *NumBuffParts = NumPartsInChunk;
  *buffPx = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffPx) != NULL);
  *buffPy = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffPy) != NULL);
  *buffPz = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffPz) != NULL);
  *buffVx = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffVx) != NULL);
  *buffVy = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffVy) != NULL);
  *buffVz = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffVz) != NULL);
  
  //read positions
  fsetpos(fp,&pos_p);
  fread(partChunk,(size_t) NumPartsInChunk,3*sizeof(float),fp);
  fgetpos(fp,&pos_p);
  for(k=0;k<NumPartsInChunk;++k)
    {
      (*buffPx)[k] = partChunk[k*3 + 0];
      (*buffPy)[k] = partChunk[k*3 + 1];
      (*buffPz)[k] = partChunk[k*3 + 2];
    }
    
  fsetpos(fp,&pos_v);
  fread(partChunk,(size_t) NumPartsInChunk,3*sizeof(float),fp);
  fgetpos(fp,&pos_v);
  for(k=0;k<NumPartsInChunk;++k)
    {
      (*buffVx)[k] = (float) (partChunk[k*3 + 0]/sqrta);
      (*buffVy)[k] = (float) (partChunk[k*3 + 1]/sqrta);
      (*buffVz)[k] = (float) (partChunk[k*3 + 2]/sqrta);
      ++CurrNumParts;
    }
  
  //clean up
  // if on last chunk, free internal buffer and set internal page num to -1 so that will init again - return 1 to indicate read of file is done
  // if not on lst call, increment chunk # and return 0 to indicate that we need to keep on reading file
  if(chunk == NumChunks-1)
    {
      assert(CurrNumParts == TotNumPageParts);
      fclose(fp);
      free(partChunk);
      InternalPageNum = -1;
      
      return 1;
    }
  else
    {
      ++chunk;
      return 0;
    }

#undef SKIP  
}

long read_GADGET_NumPages(void)
{
  FILE *fp;
  long files;
  char basefname[MAX_FILENAME];
  char fname[MAX_FILENAME];

  struct io_header_1
   {
     unsigned int      npart[6];
     double   mass[6];
     double   time;
     double   redshift;
     int      flag_sfr;
     int      flag_feedback;
     unsigned int      npartTotal[6];
     int      flag_cooling;
     int      num_files;
     double   BoxSize;
     double   Omega0;
     double   OmegaLambda;
     double   HubbleParam;
     char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
   } header1;
  int dummy;
  
#define SKIP fread(&dummy,sizeof(dummy),(size_t) 1,fp);
   
   fp = fopen(haloFinderData.SnapshotFileList,"r");
   assert(fp != NULL);
   fgets(basefname,MAX_FILENAME,fp);
   fclose(fp);
   
#ifdef DEBUG
#if DEBUG_LEVEL > 1
   fprintf(stderr,"%d: reading num parts from '%s'\n",ThisTask,basefname);
#endif
#endif
   fp = fopen(basefname,"r");
   if(fp == NULL)
     {
       sprintf(fname,"%s.%d",basefname,0);
       fp = fopen(fname,"r");
       assert(fp != NULL);
     }
   
   SKIP;
   fread(&header1,sizeof(header1),(size_t) 1,fp);
   SKIP;
  
   files = (long) (header1.num_files);
      
   fclose(fp);
#undef SKIP  

   return files;
}

